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Abstract — Continuum  kinetic  models,  such  as  Maxwell- 
Boltzmann,  present  a  viable  alternative  to  particle-in-cell  (PIC) 
models  because  they  can  be  cast  in  conservation  form  and 
are  not  susceptible  to  noise.  By  treating  the  associated  phase 
space  distribution  function  as  a  continuous  incompressible  fluid 
occupying  a  volume  of  position-velocity  space,  evolution  of  the 
distribution  function  is  determined  by  solving  a  6-D  advection 
equation.  In  cases  where  collision  terms  are  negligible,  the 
Boltzmann  model  is  reduced  to  a  Vlasov  model.  A  4th-order 
accurate  continuum  kinetic  Vlasov  model  has  been  developed 
in  one  spatial  and  one  velocity  dimension  to  address  the  chal¬ 
lenges  associated  with  solving  a  hyperbolic  governing  equation 
in  multidimensional  phase  space.  The  governing  equation  is 
cast  in  conservation  law  form  and  solved  with  a  finite  volume 
representation.  Adaptive  mesh  refinement  (AMR)  is  used  to  allow 
for  efficient  use  of  computational  resources  while  maintaining 
desired  levels  of  resolution.  Consequently,  with  AMR  the  model  is 
able  to  capture  the  fine  structures  that  develop  in  the  distribution 
function  as  it  evolves  in  time,  while  using  low  resolution  in  the 
tail  of  the  distribution  function.  The  model  is  tested  on  several 
plasma  phenomena  including:  weak  and  strong  Landau  damping 
and  the  two-stream  instability.  Conservation  properties  of  the 
method  are  Investigated. 

I.  Introduction 

In  statistical  plasma  kinetic  theory,  each  particle  species  is 
treated  as  a  distribution  function  evolving  in  position-velocity 
phase  space.  Lor  a  collisionless  plasma,  the  evolution  is 
described  by  a  coupled  system  of  partial  differential  equations: 
the  Vlasov  equation  and  Maxwell’s  equations.  When  solved  in 
its  most  general  form,  this  model  has  six  dimensions  —  three 
spatial  coordinates  and  three  velocity  coordinates  —  thereby 
making  it  computationally  costly.  Lor  this  reason,  particle-in- 
cell  (PIC)  methods,  which  rely  on  sampling  the  distribution 
function,  have  been  the  dominant  means  of  solving  this  system. 
Due  to  recent  advancements  in  supercomputing  technology, 
however,  full  phase-space  continuum  methods  have  garnered 
more  attention  and  have  become  a  viable  means  of  simulating 
nonlinear  plasma  kinetics.  [2],  [3],  [4],  [1],  [7] 

Continuum  methods  are  advantageous  because  they  can  be 
cast  in  conservation-law  form  and  thereby  be  solved  using 
flux-based  techniques  that  have  been  well-developed  in  the 
field  of  hyperbolic  partial  differential  equations.  Hyperbolic 
solvers  also  allow  for:  straightforward  parallelization  based  on 
domain  decomposition;  the  use  of  adaptive  mesh  refinement 
(AMR)  techniques;  and  numerical  methods  that  are  high-order 


accurate  in  space  and  time.  [2]  Parallel  AMR  can,  in  particular, 
be  used  to  reduce  the  cost  of  a  continuum  multi-dimensional 
Vlasov  simulation  by  focusing  computational  resources  in 
dynamic  parts  of  the  domain.  Moreover,  continuum  methods 
do  not  suffer  from  sampling-associated  noise  seen  in  PIC 
methods  [9]  and  can  thus  produce  more  physically  accurate 
results. 

The  content  of  this  paper  investigates  a  high-order  accu¬ 
rate  numerical  method  for  modeling  the  electrostatic  regime 
represented  by  the  Vlasov-Poisson  system  in  one  spatial  and 
one  velocity  dimension.  The  unsplit  finite  volume  method 
is  fourth-order  accurate  in  time  and  space  and  solves  the 
system  of  equations  in  conservation-law  form.  The  method 
is  benchmarked  against  analytic  results  for  weak  Landau 
damping,  strong  Landau  damping,  and  two-stream  instability 
evolution.  Its  conservation  properties  are  assessed  and  results 
with  adaptive  mesh  refinement  are  presented. 


II.  Vlasov-Poisson  System  and  Underlying 
Assumptions 


Plasma  evolution  timescales  over  which  electrons  are  dy¬ 
namic  and  ions  remain  static  are  considered.  In  such  a  system 
there  is  only  one  evolving  distribution  function  f{x,v,t)  — 
that  of  the  electrons.  Assuming  collisions  have  negligible 
effect  means  the  plasma  kinetics  can  be  modeled  by  the  Vlasov 
equation  [8]  in  conservation-law  form: 


m  ov 


(1) 


where  /  is  the  distribution  function,  t  is  time,  x  is  the  phase 
space  position  coordinate,  v  is  the  phase  space  x-direction 
velocity  coordinate,  e  is  the  absolute  value  of  electron  charge, 
m  is  the  electron  mass,  and  E  is  the  electric  field  in  the  x- 
direction.  In  Eq.  1  only  electrostatic  forces  are  considered  such 
that  currents  and  magnetic  fields  are  assumed  to  be  negligible. 
The  electric  potential,  (f),  is  calculated  using  Poisson’s  equation 
in  one  dimension. 


_  _Pc 
dx"^  £o 


(2) 


where  pc  is  the  charge  density  and  £o  is  the  vacuum  permit¬ 
tivity.  Note  that  the  electric  field  is  related  to  the  gradient  of 
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Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
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the  potential:  E  =  —  V(/).  The  charge  density  is  defined  as 


Pc{x)  =  e 


f{x,v,t)dv 


(3) 


Because  the  system  is  assumed  to  be  non-relativistic,  the 
range  over  which  velocity  is  defined  is  unbounded,  meaning 
in  theory  that  v  G  (—00,00).  In  practice,  i.e.  in  numerical 
simulations,  the  velocity  is  bounded  by  some  sufficiently  high 
value  VjYiax  such  that  V  G  [  The  value  of  Vxnax 

is  set  such  that  f(x,  Evmax)\t=o  machine  precision. 

The  distribution  function  f{x,  v,  t)  is  normalized  such  that  the 
net  charge  density  over  the  entire  spatial  domain  is  zero. 

The  kinetic  energy  and  potential  energy  for  the  Vlasov- 
Poisson  system  are  defined  as  follows  [5] 

KE  =^m  J  j  v'^fdvdx  (4) 

PE  =^eo  j  E'^dx.  (5) 


The  values  of  kinetic  and  potential  energy  provide  a  useful 
means  of  comparing  results  of  numerical  simulations  to  theo¬ 
retical  predictions.  The  values  of  energy  also  provide  a  metric 
by  which  to  evaluate  conservation  properties  of  the  Vlasov- 
Poisson  solver. 


III.  Discretization  and  Integration  in  Time 


A  fourth-order  finite  volume  method  is  employed  to  advance 
the  solution  f{x,v)  in  time.  This  is  done  by  first  initializing 
the  distribution  function  by  a  fourth-order  cell-average  of  the 
analytic  initial  condition,  /o,  at  discrete  set  of  points  i,j, 
denoting  the  spatial  and  velocity  index,  respectively.  Thus  at 
time  t  =  0  the  cell-average  distribution  function  to  fourth  order 
is 


(/)i.i  =foiihx,jhy)  -f 


hld^fo  ,  hld^h 


Note  that 


24  dx"^  24 
denotes  the  average  value  over  a  cell 


ifkj  = 


1 


h'v  h^] 


(6) 


f{x,v)dvdx  (7) 


and  hx  and  hy  refer  to  cell  widths  in  the  x  and  v  directions. 
The  second  derivative  terms  on  the  right  hand  side  of  Eq.  6 
are  defined  using  a  second-order  accurate  centered  difference 
stencil 

9^  fo  ~  2/ij  + 

dv^  ~  hi  ^  ’ 

9^  fo  _  fi+i,j  ~  2/i,j  + 

dx^  ~  hi  ■  ^  ’ 


The  value  of  the  cell-average  potential  {<j>)i  is  computed  from  a 
fourth-order  finite  volume  stencil  for  the  Poisson  equation  [10] 


Eq  12  d  2 

4  1 


(10) 


where  (pc)i  is  defined  in  terms  of  the  zeroth  moment  of  the 
cell-average  distribution, 

{Pe)i  =  l-hyY,{f)^,J■  (11) 

j 

The  cell-average  electric  field  is  computed  from  the  potential: 

To  advance  the  solution  {f)ij  in  time,  it  is  necessary  to 
calculate  the  distribution  function  flux  entering  each  cell  face 
of  cell  i,j.  Based  on  the  velocity  and  the  electric  field,  the 
fourth-order  accurate  fluxes  are  calculated  as  follows: 


+  (13) 

=  (/)., +  i 

+  ^{k)^+ld+^  -  (14) 


Once  the  values  of  flux  at  each  cell  face  are  known,  the  Vlasov 
equation  is  reduced  to  an  ordinary  differential  equation: 


dt 


1 

hx 

e  1 

m  h, 


(15) 


Eq.  15  is  integrated  in  time  using  a  fourth-order  explicit 
Runge-Kutta  algorithm  to  advance  {f)ij  by  a  time  step. 
The  updated  cell-average  distribution  function  is  then  used  to 
compute  the  electric  field,  and  thus  to  advance  the  solution. 
This  amounts  to  a  fourth-order  accurate  discretization  in  time 
and  space,  using  centered  differencing  for  the  latter. 


IV.  Single-Grid  Simulation  Results 


The  single  grid  algorithm  described  in  the  previous  section 
is  benchmarked  against  analytic  results  from  kinetic  theory. 
A  standard  test  case  is  to  simulate  weak  Landau  damping  by 
initializing  a  Maxwellian  distribution  in  velocity  space  with  a 
small  position-dependent  perturbation: 


exp 


(1  -I-  acos{kx)) , 


(16) 


with  k  =  0.5  and  a  =  0.01.  In  weak  Landau  damping, 
potential  energy  (see  Eq.  5)  is  transformed  into  thermal  energy 
as  indicated  by  a  steady  net  decrease  in  the  value  of  the  former. 
The  simulated  decay  rate  of  potential  energy  in  time  is  in 
good  agreement  with  the  theoretical  prediction,  as  shown  in 
Fig.  1(a). 

To  test  a  non-linear  regime  of  plasma  kinetics,  the  simula¬ 
tion  is  benchmarked  against  theoretical  predictions  for  strong 
Landau  damping.  In  strong  Landau  damping,  the  Maxwellian 
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Normalized  Time  (t  =103) 

'  n  p' 

(a)  Weak  Landau  damping 


Fig.  1;  Simulation  results  for  (a)  weak  Landau  damping  and 
(b)  strong  Landau  damping  using  a  256  x  256  grid.  The 
evolution  of  potential  energy  demonstrates  close  agreement 
with  theoretical  predictions. 

in  Eq.  16  is  given  a  large  perturbation  with  amplitude  a  =  0.5. 
In  this  case,  the  potential  energy  in  the  system  evolves  non- 
linearly.  This  evolution  is  shown  in  Fig.  1(b),  where  the  decay 
and  growth  rates  in  potential  energy  from  the  simulation  are 
shown  in  comparison  with  theoretical  predictions. 

The  algorithm  is  also  used  to  simulate  the  two-stream 
instability.  This  simulation  is  initialized  using  the  following 
distribution  function 

/(a;,i^)lt=o  =:^^^exp  (l  +  acos(/ca;))  (17) 

with  k  =  0.5  and  a  =  0.01.  The  spatial  perturbation  results  in 
inhomogeneities  in  the  electron  distribution  such  that  kinetic 
energy  is  converted  into  potential  energy.  Fig.  2  shows  how 


Fig.  2:  Evolution  of  the  two-stream  instability  for  normalized 
time  tn  =  0, 18.8,  20.4,  29.8.  Note  that  at  approximately  t  = 
18.8  the  instability  becomes  non-linear. 


Fig.  3:  The  two-stream  instability  results  demonstrate  close 
agreement  with  theoretical  predictions  for  the  instability 
growth  rate. 


the  electrons  become  trapped,  as  is  evidenced  by  the  formation 
of  a  vortex-like  structure  in  the  phase  space  distribution.  It  is 
important  to  note  that  dispersive  effects,  which  can  be  seen 
at  time  t  —  29.8,  cause  the  distribution  function  to  become 
negative  in  certain  parts  of  the  domain.  This  is  an  unphysical 
result  that  can  be  resolved  through  the  use  of  limiters.  When 
compared  to  linear  theory,  the  growth  rate  in  the  potential 
energy  resulting  from  the  simulation  is  consistent  with  the 
theoretically  predicted  value  for  the  two-stream  instability,  as 
shown  in  Fig.  3. 
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Fig.  4:  Two-stream  instability  at  t  —  22.1.  AMR  simulation 
with  64  X  64  base  grid,  three  levels  of  grid  rehnement,  and 
a  rehnement  ratio  of  four.  Thus  the  smallest  boxes  indicate 
regions  of  the  domain  that  have  sixteen  times  the  resolution 
of  the  coarsest  grid. 


V.  Two-stream  Instability  with  AMR 

Adaptive  mesh  rehnement  is  implemented  into  the  single¬ 
grid  algorithm  described  above,  using  the  techniques  described 
in  Ref.  [6].  As  indicated  in  Eq.  3,  the  charge  density  is 
computed  by  taking  a  velocity  moment  of  the  distribution 
function  at  each  spatial  location.  To  do  this  on  multiple  levels 
of  grids  requires  that  additional  interpolation  steps  be  added 
to  the  AMR  algorithm. 

First  a  coarse-grid  number  density  is  computed  using  the 
coarse  cell-average  distribution  function.  The  coarse  grid 
number  density  is  interpolated  to  the  hne  grid.  In  order  to 
get  the  most  accurate  value  of  number  density,  the  integral 
contribution  from  hne  grid  must  be  included.  This  is  done  by 
taking  the  difference  between  the  coarse  and  hne  distribution 
functions  in  regions  where  the  hne  grid  exists,  and  integrating 
this  difference  with  respect  to  v.  This  position-dependent 
contribution  is  then  added  to  the  interpolated  coarse  grid 
number  density  in  spatial  regions  where  the  hne  grid  exists. 

As  a  demonstration  of  adaptive  mesh  rehnement  capability, 
the  two-stream  instability  is  simulated  using  a  three-level 
adaptively  rehned  grid.  The  distribution  function  and  grid  are 
shown  at  time  f  =  22.1  in  Fig.  4.  The  simulation  uses  a  base 
grid  of  size  64  x  64  with  a  factor  of  four  rehnement  for  each 
successive  level.  Further  analysis  is  required  to  evaluate  the 
merits  of  AMR  and  to  determine  to  what  extent  it  reduces  the 
computational  cost  of  Vlasov-Poisson  simulations. 
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(c)  Energy  change  relative  to  total  energy  at  t  =  0 

Fig.  5;  Absolute  values  of  mass,  momentum  and  energy  for 
weak  Fandau  damping.  Mass  is  conserved  to  10“^^  precision, 
momentum  is  conserved  to  machine  precision,  and  total  energy 
conservation  depends  on  the  grid  resolution. 
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VI.  Conservation  Properties 


References 


The  conservation  properties  of  the  described  algorithm  are 
investigated  by  tracking  the  value  of  the  zeroth,  first,  and 
second  velocity  moment  of  the  distribution  function  (inte¬ 
grated  over  the  entire  domain)  as  a  function  of  time.  These 
integrals  amount  to  conservation  of  mass,  momentum  and 
energy.  Theoretically  the  value  of  each  should  be  conserved 
to  machine  precision,  but  in  practice  numerical  errors  are 
introduced  through  the  fourth-order  approximation  to  the 
partial  differential  equation.  Fig.  5  shows  to  what  extent  the 
mass,  momentum,  and  energy  of  the  system  are  conserved 
for  the  case  of  weak  Landau  damping.  The  zeroth  moment 
is  preserved  to  10“^^  independent  of  the  grid  resolution.  The 
first  moment  is  zero  to  machine  precision,  and  its  value  is  also 
independent  of  grid  resolution.  The  second  moment,  as  mea¬ 
sured  by  the  change  in  total  energy  relative  to  the  initial  value 
at  f  =  0,  attains  a  value  that  is  several  orders  of  magnitude 
greater  than  machine  precision.  The  second  moment  remains 
bounded  for  all  time,  and  fourth-order  convergence  of  the 
numerical  method  is  demonstrated  by  the  fact  that  the  change 
in  total  energy  decreases  by  a  factor  of  approximately  sixteen 
as  the  grid  resolution  is  doubled  in  x  and  v.  Conservation 
properties  for  non-linear  plasma  phenomena  require  further 
investigation. 


[1]  T.  D.  Arber  and  R.  G.  L.  Vann.  “A  critical  comparison  of  Eulerian-grid- 
based  Vlasov  solvers,”  Journal  of  Computational  Physics,  Vol.  180,  No. 
1,  pp.  339-357,  June  2002. 

[2]  1.  W.  Banks,  and  1.  A.  F.  Hittinger.  “A  new  class  of  nonlinear  finite- 
volume  methods  for  Vlasov  simulation,”  IEEE  Transactions  on  Plasma 
Science,  Vol.  38,  No.  9,  September  2010. 

[3]  N.  V.  Elkina  and  1.  Buchner.  “A  new  conservative  unsplit  method  for  the 
solution  of  the  Vlasov  equation,”  Journal  of  Computational  Physics,  Vol. 
213,  No.  2,  pp  862-875,  April  2006. 

[4]  F.  Filbet,  E.  Sonnendrucker.  “Comparison  of  Eulerian  Vlasov  solvers,” 
Computer  Physics  Communications,  Vol.  150,  No.  3,  pp.  247-266,  June 
2002. 

[5]  R.  T.  Glassey.  Cauchy  Problem  in  Kinetic  Theory,  Society  of  Industrial 
and  Applied  Mathematics,  Philadelphia,  1996. 

[6]  P.  McCorquodale  and  P.  Colella.  “A  High-Order  Finite- Volume  Method 
for  Conservation  Laws  on  Locally  Refined  Grids,”  Communications  in 
Applied  Mathematics  and  Computational  Science,  Vol.  6,  No.  1,  pp.  1- 
25,  2011 

[7]  N.  J.  Sircombe  and  T.  D.  Arber.  “VALIS:  A  split-conservative  scheme 
for  the  relativistic  2D  Vlasov-Maxwell  system,”  Journal  of  Computational 
Physics,  Vol.  228,  No.  13,  pp  4773-4788,  July  2009 

[8]  N.  A.  Krall  and  A.  W.  Trivelpiece.  Principles  of  plasma  physics,  McGraw- 
Hill,  1973. 

[9]  B.  Wang,  G.  H.  Miller,  and  P.  Colella.  “A  Particle-in-cell  method  with 
adaptive  phase-space  remapping  for  kinetic  plasmas,”  SIAM  Journal  of 
Scientific  Computing,  Vol.  33,  No.  6,  pp.  3509-3537,  September  2011. 

[10]  Q.  Zhang,  H.  Johansen,  and  P.  Colella.  “A  fourth-order  accurate  finite- 
volume  method  with  structured  adaptive  mesh  refinement  for  solving  the 
advection-diffusion  equation.”  SIAM  Journal  of  Scientific  Computing,  Vol. 
34,  No.  2,  pp.  B179-B201,  2012 


VII.  Conclusions 

A  fourth  order  accurate  algorithm  in  space  and  time  has 
been  developed  to  solve  the  Vlasov-Poisson  system  in  one 
spatial  and  one  velocity  dimension.  The  simulation  results 
demonstrate  close  agreement  with  theoretical  predictions,  as 
tested  on  weak  Landau  damping,  strong  Landau  damping,  and 
the  two-stream  instability.  AMR  is  also  demonstrated  as  a 
practical  means  of  reducing  computational  cost  by  focusing 
computational  resources  in  dynamic  regions  of  the  domain. 
Further  work  needs  to  be  done  to  assess  the  speed-up  offered 
by  adaptive  mesh  refinement,  specifically  as  it  applies  in  multi¬ 
ple  dimensions.  The  described  algorithm  will  also  benefit  from 
the  use  of  limiters,  which  will  address  dispersive  oscillations 
and  lack  of  positivity  preservation  that  are  inherent  to  high- 
order  finite  volume  calculations. 
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